Influence of pore-scale disorder on viscous fingering during drainage 
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We study viscous fingering during drainage experiments in linear Hele-Shaw cells filled with a 
random porous medium. The central zone of the cell is found to be statistically more occupied 
than the average, and to have a lateral width of 40% of the system width, irrespectively of the 
capillary number Ca. A crossover length Wf oc Ca~^ separates lower scales where the invader's 
fractal dimension D ~ 1.83 is identical to capillary fingering, and larger scales where the dimension 
is found to be D ~ 1.53. The lateral width and the large scale dimension are lower than the results 
for Diffusion Limited Aggregation, but can be explained in terms of Dielectric Breakdown Model. 
Indeed, we show that when averaging over the quenched disorder in capillary thresholds, an effective 
law i; oc (VP)^ relates the average interface growth rate and the local pressure gradient. 

PACS numbers: 47.20.Gv, 47.53.+n, 47.54.+r, 47.55.-t, 47.55.Mh, 68.05.-n, 68.05.Cf, 81.05.Rm. 



Viscous fingering instabilities in immiscible two-fluid 
flows in porous materials have been intensely studied 
over the past 50 years 0, both because of their impor- 
tant role in oil recovery processes, and as a paradigm of 
simple pattern forming system. Their dynamics is con- 
trolled by the interplay between viscous, capillary and 
gravity forces. The ratio of viscous forces to the cap- 
illary ones at pore scale is quantified by the capillary 
number Ca — /iw/a^/(7«;), where a is the characteristic 
pore size, vf is the filtration velocity, 7 the interfacial 
tension, k the permeability of the cell, and the viscos- 
ity of the displaced fluid, supposed here much larger than 
the viscosity of the invading one. 

There is a strong analogy between viscous fingering in 
porous media and Diffusion Limited Aggregation (DLA) , 
as was first pointed out by Patcrson 2]. Indeed, both pro- 
cesses of DLA and viscous fingering in empty Hclc-Shaw 




FIG. 1: Invasion clusters on thresholded images at capil- 
lary numbers Ca = 0.06 (a) and 0.22 (b,c), for W/a = 210 
(a,b)and 110 (c), with displayed system lateral boundaries. 
The superimposed gray map shows the occupancy probability 
Tv{x, z) of the invader, in a moving reference frame attached to 
the most advanced invasion tip and to the lateral boundaries. 



cells belong to the family of Laplacian Growth Models, 
i.e. obey the Laplacian growth equation V^P = 0, with 
an interfacial growth rate v oc —VP, where P is the dif- 
fusing field, i.e. the probability density of random walk- 
ers in DLA, or the pressure in viscous fingering. Despite 
differences as respectively a stochastic and deterministic 
growth, and boundary conditions as respectively P — Q 
or P = — 7/r with r the interfacial curvature, it is admit- 
ted that these processes belong to the same universality 
class 1^ H, In radial geometry these processes lead 
to fractal structures of dimensions D = 1.70 ± 0.03 |^, 
1.713 ± 0.0003 0, and 1.7 respectively in viscous fin- 
gering in empty Hele-Shaw cells, DLA, and numerical 
solutions of deterministic Laplacian growth. The two 
numerical models have been reexplored recently using 
stochastic conformal mapping theory However, 
in Hclc-Shaw cells filled with disordered porous materi- 
als similar to the one used here, a lower fractal dimension 
D = 1.58 ± 0.09 has been measured 0- 

In straight channels, DLA gives rise to fractal struc- 
tures of dimension 1.71, occupying on average a lateral 
fraction A = 0.62 of the system width W 8]. Viscous 
fingering in empty Hele-Shaw cell converge towards the 
Saffman Taylor (ST) solution 1^, with a uniformly prop- 
agating fingerlike interface covering a fraction A = 0.5 
of the system width at large capillary numbers 0, 0| , 
selected by the interfacial tension ■ 

In the system we study, the cell is filled with a disor- 
dered porous medium, and the non- wetting invader of low 
viscosity shows a branched structure that depends on Ca 
(Fig.^. We will show theoretically that if indeed at high 
capillary number, the process is well described by DLA 
as often suggested 0, there is also at intermediate Ca a 
regime where the flow in random porous media is better 
described by another Laplacian model, namely a Dielec- 
tric Breakdown Model (DBM) with rj = 2~ the interfacial 
growth rate is w oc (VP)'' in DBM, rj — I corresponding 
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FIG. 2: Scaled invader's density n{z)/nca as function of 
2z/W , distance to the most advanced tip scaled by the sys- 
tem's halfwidth, for three system sizes and four capillary num- 
bers. The lines are the experimentally determined cumulative 
growth probability, and the theoretical Saffman Taylor solution 
for a single finger occupying a lateral fraction lambda = 0.4 
of the system. 



to DLA -. Our conclusion is supported by both experi- 
mental evidence based on a characterization of the large- 
and small-scale geometry of the invader, and theoretical 
arguments based on averaging the capillary forces con- 
tribution to the interface growth rate over the pore scale 
randomness. The vanishing capillary number limit of our 
system corresponds to capillary fingering, where pores 
are invaded one by one, forming a propagating front fill- 
ing the whole system, leaving behind isolated clusters of 
viscous fluid. The invader's fractal dimension is D = 1.83 
[l^ . theoretically explained by the invasion percolation 
model At moderate but finite capillary numbers un- 
der interest here, we will show that small scales still cor- 
respond to invasion percolation, up to a crossover scale 
Wf/a oc Ca~^ characterizing the maximum size of iso- 
lated clusters of defending fluid, above which the system 
geometry can be described by the DBM with r/ = 2. 

We study viscous fingering processes in linear Hele- 
Shaw cells of thickness a = 1mm filled at 38% with 
a monolayer of randomly located immobile glass beads 
of diameter a, between which air displaces a solution 
of 90% glycerol - 10% water of much larger viscosity 
/i — 0.165 Pa-s, wetting the beads and walls of the cell, 
i.e. in drainage conditions. The interfacial tension and 
the permeability of the cell are respectively 7 = 0.064 
N-m~i and k = 0.00166 ± 0.00017 mm^. We investi- 
gate regimes ranging from capillary to viscous fingering 
(0.01 < Ca < 0.5), in cells with impermeable lateral 
walls and dimensions W x L x a, with widths perpen- 
dicularly to the fiow direction W/a = 110, 215 and 430, 
and a length L/a — 840. The cell is set horizontally, so 
that gravity is irrelevant. A constant filtration rate of 
water-glycerol is ensured by a controlled gravity-driven 
pump. 

Pictures of the flow pattern arc taken from the top, and 
treated to extract the invading air cluster (with pixels of 
size 0.55a), as the black clusters in Fig.^ In rcf. 14], we 



have shown that the invasion process is stationary, up to 
fluctuations arising from the disorder in pore geometries. 
To extract the underlying average stationary behavior, 
all quantities are then analyzed in the reference frame 
{x, z) attached to the lateral boundaries at x — and 
X = W, and to the foremost propagating tip at z = 0, z 
pointing against the flow direction (this tip indeed prop- 
agates at a roughly constant speed vup llJ|)- Average 
quantities at any position {x, z) of the tip related frame, 
are defined using all stages and points of the invasion pro- 
cess, excluding regions closer than W/2 from the inlet or 
outlet, to avoid finite size effects. 

The average occupancy map Tr{x, z) is defined as [l4| : 
for each time (or each picture), we assign the value 1 to 
the coordinate (x, z) if air is present there, otherwise. 
7r(x, z) obtained as the time average of such occupancy 
function, is displayed as graymap in Fig. ^ 

Next we compute the average number of occupied 
pores per unit length at a distance z behind the tip, n{z)^ 

which is related to tt as n{z) = (1/a^) /J^ 7r(a;, z) dx. We 
show in Fig.Ela data collapse for different capillary num- 
bers Ca and system widths W, n{z)/nca ~ $(2z/VF), 
where nca = {W/a^)vf /vup 0|. $ is a function increas- 
ing from at z = towards 1 at z = -l-oo, as granted 
by conservation of the displaced fluid for a statistically 
stationary process 0|. <I> evaluated in Fig. El is a cu- 
mulative growth probability defined in Ref. jl4| . and is 
obtained as an average over all experiments and sizes. 

We also characterize the lateral structure of the in- 
vader in the frozen zone, z > where less than 
10% of the invasion activity takes place since $(2) > 
0.9. We define over this zone a distribution p{x) = 
\W/(a'^nca)] 7r(a;, 00) — (t'tip/wf) 7r(a;,oo), so that 
/J^ p{x)dx/W ~ 1. This quantity, presented in Fig.lSfa) 
for an average over five experiments at capillary numbers 
Ca = 0.06 and 0.22 for W/a = 215, and over four exper- 
iments with 0.06 < Ca < 0.22 for W/a — 110, is reason- 
ably independent of the capillary number and the system 
size, though the noise is larger at highest speed. The frac- 
tion A of the system, occupied by the invader at satura- 
tion is evaluated as in 0: A = i/pmax, or alternatively 
A = (a;+ - x~)/W, where p{x+) = p{x~) = Pmax/'^- 
Both definitions lead to A ~ 0.4 ± 0.02 for the capillary 
numbers probed, as shown in Fig. |3fa), which is signif- 
icantly smaller than the off-lattice DLA result A = 0.62 

The 2D occupancy map ^{x, z) itself, displayed as 
graymap in Fig. [T] has a maximum TTmax = Pmaxi'//'ftip, 
along a line at {x = W/2, z > W). Similarly to Arneodo's 
procedure 0, we determine the support of tt > TTmax/2, 
displayed in Fig.|3{b) and (c), which corresponds to the 
most often occupied region. Within the noise error, the 
shape of this region resembles the theoretical ST curve 
corresponding to this A = 0.4 (gray lines in Fig. 
For such ST finger, this curve would also correspond to 
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FIG. 3: (a) Normalized occupation density p{x), with half- 
maximum reached over a width QAW . (b,c) Average occupa- 
tion density map of the invader thresholded at half-maximum, 
for a system size W/a—215, in the reference frame attached 
to the tip position, at Ca = 0.06 (b) and 0.22 (c), compared 
to ST curves for X = 0.35 and 0.45. 
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FIG. 4: (a) Mass fractal dimension of the invasion cluster at 
various capillary numbers, (b): Cross-over scale wf between 
capillary and viscous scales, as function of Ca. 



the scaled longitudinal density of invader ^(z). $(2) de- 
termined for our experiments (continuous line in Fig. \^ 
and this theoretical ST shape (dashed line) also present 
some similarities. Systematic deviations from the math- 
ematical ST solution at A = 0.4 might nonetheless exist, 
since they have been seen between the DLA envelopes 
and the corresponding ST solution at A = 0.62 [l5| . 

The mass fractal dimension of the invasion clusters is 
analysed by box-counting. N{s) is the number of boxes 
of size s to cover the invader. Fig. ^displays a normalized 
distribution N{s)/N{a/Ca) as function of (s/a) • Ca for 
various capillary numbers. By linear regression of this 
collapsed log-log data, we find that N{s) ^ g-i-83±o.oi 
for small scales s < a/Ca, and N(s) ~ g-i-53±o.o2 j-^j, 
larger scales. The result can be explained by the fol- 
lowing approximations. The distribution of pore throats 
sizes results in a distribution of capillary pressure thresh- 
olds Pt, g{Pt), of characteristic width Wt- Consider a box 
of size Wf along the cluster boundary in the active zone, 
such that Wf ■ VPf, = Wt where VPb is a characteris- 
tic pressure gradient. At scales s < Wf, viscous pressure 



variations are lower than capillary threshold fluctuations, 
and the most likely invaded pores correspond to the low- 
est random thresholds, which corresponds to capillary 
fingering, thus leading to D = 1.83. Conversely, at larger 
sizes (s > Wf), the invasion activity is determined essen- 
tially by the spatial variations of VPf,. Assuming that 
VPfc scales as the imposed VP(— 00) oc Ca, i.e. neglect- 
ing the geometry variations between different speeds, Wf 
scales as Wt/VFf, oc a/Ca, as confirmed by the data 
collapse in Fig.^fa). Wf can also be determined experi- 
mentally as a characteristic branch width, since capillary 
fingering leaves isolated clusters of trapped fluid, whereas 
the large scale structure is branched. After removing all 
trapped clusters, we determine Wf as the average length 
of intersects of the structure from cuts along x. Indeed, 
Fig. ^b) is consistent with a scaling law Wf /a oc Ca~^, 
below a saturation at large Ca. 

Eventually, we sketch a possible explanation for the 
width selection A — 0.4 and the large scale fractal di- 
mension D ~ 1.53 ± 0.02, which are smaller than their 
counterparts in DLA, respectively 0.62 and 1.71. Ne- 
glecting the small scale permeability variations leads to 
a Laplacian pressure field in the defending fluid. The 
boundary condition for the pressure field is then \/P(z = 
—00) — —fivf/n, and VP(a; = 0,14^) ■ x = where x 
is the unit vector along x. The dynamics of the pro- 
cess is then entirely controlled by the boundary condi- 
tion along the invading fluid, i.e. by the capillary pres- 
sure drop across the meniscus in the pore neck and the 
pressure gradient in the invaded fluid. For a given pres- 
sure difference at pore scale between the invading air at 
Pq, and the pressure Pi in the glycerol- filled pore, we 
decompose Pq — Pi = APy -\- Pc, where APy is a vis- 
cous pressure drop in the pore neck, and the capillary 
pressure drop is Pc = j/r -\- 2^/a, where the in- and 
out-of-plane curvature of the interface are respectively r 
and a/ 2. As a meniscus progresses between neighbor- 
ing beads, its curvature goes through a minimum in 
the pore neck. The meniscus will be able to pass the 
neck if the pressure drop Pq — Pi exceeds the thresh- 
old Pt — Pc{rm)- For the sake of simplicity, the prob- 
ability distribution of the thresholds g{Pt) is considered 

flat, between Pmin and Pmax, with Wt = Pmax - Pmin 

and giPt) = e{Pt - Prmn)0{Pmax ^ Pt)lWt, whcrc 6 is 
the Heaviside function. In the pure capillary fingering 
limit Ca 0, the pressure field P is homogeneous in 
the defending fluid, and a pore is invaded when Pq — P 
reaches the minimum threshold along the boundary, close 
to Pmin- At higher capillary number, we want to relate 
the invasion rate to the local capillary threshold, and to 
the pressure Pi in the liquid-filled pore nearest to the in- 
terface. If Pq — Pi < Pt, the meniscus adjusts reversibly 
in the pore neck, and the next pore is not invaded. Con- 
versely, if Po ^ Pi > Pt, the pore will be invaded, and 
most of the invasion time is spent in the thinnest region 
of the pore neck. A characteristic interface velocity can 
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be evaluated by the Washburn equation 16] at this point: 

V {2K/na){Po -Pi- Pt)0{Po - Pi - Ft), where the 

heaviside function results from a zero invasion velocity 
if the pore is not invaded. Hypothesizing that only the 
average growth rate controls the process, independently 
of the particular realization of random thresholds, the 
growth rate averaged over all possible pore neck config- 
urations, is 

< >= / — (Po -Pi- Pt)0{Po -Pi- Pt)9(Pt)dPt 
= -{n/tia)0{Po - Pi - P„„„) X (1) 

{[(PO - Pi - P,mn)VWt]e[P^a. " (Pq " ^l)] + 
2[Po — Pi — {Pmin + Pmax)/'2']0[Po — Pi — Pmax]} ■ 

At moderate capillary numbers, such as Pq — Pi < Pmax, 
if we assume that the capillary pressure drop is around 
Pc — Pmin when the invasion meniscus is at the entrance 
of the pore neck, we note that (Pq — Pmin — Pi)/ct = 
APy/a ~ VP/2, and Eq. implies that the growth 
rate goes as < w >= — ok/(4/iM^()(VP)^. This effective 
quadratic relationship between the average growth rate 
and the local pressure gradient arises from the distribu- 
tion of capillary thresholds, and means that such inva- 
sion process should be in the universality class of DBM 
with 77 = 2, rather than DLA (DBM, 77 = 1). Indeed, 
in DBM simulations in linear channels, A is a decreas- 
ing function of rj (as in related deterministic problems, 
as viscous fingering in shear-thinning fluids |l7l |. or 77- 
model [Till), and Somfai et al. 15] report A ~ 0.62 and 
0.5 for respectively rj ~ 1 and 1.5, so that the observed 
A = 0.4 is consistent with 77 = 2. The fractal dimension 
of DBM is also a decreasing function of 77, and rj = 2 cor- 
responds D = 1.4 ±0.1 which is close to the observed 
D = 1.53 ± 0.02 in our experiments. 

Note that at high capillary numbers such that locally 
Pq— Pi ^ Pmax, thc threshold fluctuations are not felt by 
the interface, and Eq. ^ leads to < >— —{K/fi)\7P, 
which would correspond to a classic DLA process. We 
have checked by numerically solving the Laplace equa- 
tion with the experimental clusters as boundaries that 
all experiments performed here were at moderate enough 
capillary number to have Pq — Pi < Pmax all along 
the boundary |1J|, i.e. the quadratic law < v >= 
-aK/(4^Wt)(VP)2 is expected to hold. 

Even at moderate Co, deviations from the DBM model 
with 77 = 2 could be observed for significantly non-flat 
distribution of the capillary thresholds in the random 
porous medium, for which Eq.l^ would lead to a more 
complicated dependence of the growth rate v on VP, re- 
flecting the details of this distribution, and not simply a 



power-law effective relationship. It would be interesting 
in future work to explore numerically and experimen- 
tally the detailed effect of non-flat capillary threshold 
distributions on the selected fractal dimension, average 
width occupied in the system, and total displaced mass 
ncaiCa) (reported in [ij] for the present work), to ex- 
tract the influence of the disorder on the best capillary 
number to select in order to maximize the efficiency of 
the extraction process. 
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